* _______________ HOUSEKEEPING _________________________________________________

	//	Specify version
	
		version 13
	
	//	Clear everything, change working directory, and set more off

		cls
		clear
		macro drop all
		cd "/Users/davidflood/Desktop"
		graph drop _all
		set more off

	//	Use children's recode .dta file	
		
		use "/Users/davidflood/Dropbox/Guatemala/DHS and ENSMI/Analysis/DHS 2014-2015/Dataset and labels/GUKR71FL_copy.DTA"
	
	//	SSC packages used in producing this .do file
	
		*	ssc install postrcspline
		*	ssc install mdesc
		
	//	Log work
		
		log using "/Users/davidflood/Dropbox/Guatemala/DHS and ENSMI/Analysis/DHS 2014-2015/Logs/log_v28.log", replace
 
	//	Stata sampling procedures for DHS
	
		/* See: https://www.youtube.com/watch?v=YpXPWMUsb94
				general form: svyset [pw=weight], psu(psu) strata(strata)
				PSU = primary sampling unit = cluster (or "sector" in Guate DHS)
				strata = stratification (department x rural/urban in Guate DHS)	*/
				
		//	Renaming sampling variables for ease of coding and interpretation
		
			rename v021 cluster
			rename v022	strata
				
		//	Create weight variable and divide by 1,000,000
			
			rename v005 survey_weight_raw
			gen survey_weight = survey_weight_raw/1000000
		
		//	Set survey parameters
			
			svyset cluster [pw=survey_weight], strata(strata) singleunit(centered)
				
		//	Describe survey
			
			svydescribe survey_weight
							
		//	Describe dataset

			describe, short

			
* _______________ DEPENDENT VARIABLES __________________________________________
						
	
//	Adjust child height-for-age Z-score by dividing by 100
	
	gen haz = .
	replace haz = hw5/100
	replace haz = . if missing(haz) | hw5>9990
	label variable haz "Height-for-age Z-score"
	
//	Generate stunting indicators

	gen haz_stunted = .
	replace haz_stunted = 1 if haz <= -2
	replace haz_stunted = 0 if haz > -2
	replace haz_stunted = . if missing(haz)
	label variable haz_stunted "Stunted"
	label define haz_stunted_label 0 "Not stunted" 1 "Stunted"
	label values haz_stunted haz_stunted_label
	
	gen haz_sev_stunted = .
	replace haz_sev_stunted = 1 if haz <= -3
	replace haz_sev_stunted = 0 if haz > -3
	replace haz_sev_stunted = . if missing(haz)
	label variable haz_sev_stunted "Severely stunted"
	label define haz_sev_stunted_label 0 "Not severely stunted" 1 "Severely stunted"
	label values haz_sev_stunted haz_sev_stunted_label
	
	
* _______________ INDEPENDENT VARIABLES ________________________________________


//	Independent variables for models on contraceptive use and choice

	//	Generating prior contraceptive method
		
		gen method_prior = 0
		replace method_prior = 1 if inlist(v359,8,9)
		replace method_prior = 2 if inlist(v359,1,3,4,5,7,13,14,15,16,17,18,19) | inlist(v359,2,6,11)
		label variable method_prior "Contraceptive method, prior"
		label define method_prior_label 0 "None or folkloric" 1 "Traditional" 2 "Modern"
		label values method_prior  method_prior_label
		
	//	Generating current contraceptive method
		
		gen method_current = 0
		replace method_current = 1 if inlist(v312,8,9)
		replace method_current = 2 if inlist(v312,1,3,4,5,7,13,14,15,16,17,18,19) | inlist(v312,2,6,11)
		label variable method_current "Contraceptive method, current"
		label define method_current_label 0 "None or folkloric" 1 "Traditional" 2 "Modern"
		label values method_current method_current_label		

	//	Generating duration of use of current method

		rename v337 months_method_cont
		gen duration_method_cat = .
		replace duration_method_cat = 0 if months_method_cont >=0 & months_method_cont <15
		replace duration_method_cat = 1 if months_method_cont >=15
		replace duration_method_cat = 99 if missing(months_method_cont)
		label variable duration_method_cat "Duration of method"
		label define duration_method_cat_label 0 "<15 months" 1 "15+ months" 99 "No current use"
		label values duration_method_cat duration_method_cat_label			
		
	//	Generating current method with duration
		
		gen method_current_duration = .
		replace method_current_duration = 0 if method_current == 0
		replace method_current_duration = 1 if method_current == 1
		replace method_current_duration = 2 if method_current == 2 & duration_method_cat == 0
		replace method_current_duration = 3 if method_current == 2 & duration_method_cat == 1
		label variable method_current_duration "Contraceptive method, current"
		label define method_current_duration_label 0 "None or folkloric" 1 "Traditional" 2 "Modern, <15 months" 3 "Modern, 15+ months"
		label values method_current_duration method_current_duration_label
	
//	Independent variables for models on unmet need
	
	* Note: For details on unmet need, see https://dhsprogram.com/topics/unmet-need.cfm
	* Note 2: see p. 11 of Bradly/DHS (2012) for quick summary table and figure
	* Note 3: This DHS survey calculates unmet need also among unmarried women 
	
	//	Overall unmet variable
	
		/* 	Note that v626a is the corrrect variable for the revised definition
			See also: https://dhsprogram.com/topics/unmet-need.cfm	*/
	
		gen unmet = .
		replace unmet = 1 if inlist(v626a,1,2)
		replace unmet = 0 if inlist(v626a,0,3,4,5,6,7,8,9)
		label variable unmet "Unmet need binary"
		label define unmet_label 0 "No unmet need" 1 "Unmet need"
		label values unmet unmet_label
		
	//	Pregnancy intention
		
		rename m10 pregnancy_intention
		
		/*	1 "Then"
			2 "Later"
			3 "No more"	*/

		
* _______________ CONFOUNDING VARIABLES ________________________________________


//	Child sex

	rename b4 child_sex
	
//	Child age in months

	rename hw1 age_child_months
	gen age_child_month_range = .
	replace age_child_month_range = 1 if inrange(age_child_months,0,12)
	replace age_child_month_range = 2 if inrange(age_child_months,13,24)
	replace age_child_month_range = 3 if inrange(age_child_months,25,36)
	replace age_child_month_range = 4 if inrange(age_child_months,37,48)
	replace age_child_month_range = 5 if inrange(age_child_months,49,60)
	
// 	Education variables

	//	Recoding education labels
	
		label define education_cat_label ///
		0 "No education" 1 "Incomplete primary"	2 "Complete primary" 3 "Incomplete secondary" 4 "Complete secondary" 5 "Higher" 99 "Missing"
	
	//	Maternal education
	
		rename v149 maternal_education_cat
		label variable maternal_education_cat "Maternal education attainment"
		label values maternal_education_cat education_cat_label
	
	//	Paternal education
		
		rename v729 paternal_education_cat
		replace paternal_education_cat = 99 if paternal_education_cat == 8 | missing(paternal_education_cat) 
		label values paternal_education_cat education_cat_label
		
//	Maternal age
	
	rename v447a mother_age
	
//	Wealth index	

	rename v191 wealth_index_cont
		
//	Maternal literacy

	gen maternal_literacy = .
	replace maternal_literacy = 0 if inlist(v155,0,3)
	replace maternal_literacy = 1 if inlist(v155,1)
	replace maternal_literacy = 2 if inlist(v155,2)
	label variable maternal_literacy "Maternal literacy"
	label define maternal_literacy_label 0 "Not literate" 1 "Semi-literate" 2 "Literate"
	label values maternal_literacy maternal_literacy_label

//	Region of country

	rename v024 region
	
//	Area of residence

	rename v102 rural_urban
	
//	Generate indigenous ethnicity indicator, using auto-identification
		
	gen ethnicity = .
	replace ethnicity = 1 if inlist(setid,1)
	replace ethnicity = 0 if inlist(setid,0)
	label variable ethnicity "Indigenous dichotomous"
	label define ethnicity_label 0 "Not indigenous" 1 "Indigenous"
	label values ethnicity ethnicity_label

//	Generate marital status indicator
	* "partnered" defined here as "married" or "living with partner"

	gen marital_status = .
	replace marital_status = 0 if inlist(v501,0)
	replace marital_status = 1 if inlist(v501,1,2)
	replace marital_status = 2 if inlist(v501,3,4,5)
	label variable marital_status "Marital status"
	label define marital_status_label 0 "Never in union" 1 "Current partner" 2 "Formerly had partner"
	label values marital_status marital_status_label
				
//	Generate diarrhea indicators (if diarrhea in last 2 weeks)

	gen diarrhea = .
	replace diarrhea = 1 if inlist(h11,1,2)
	replace diarrhea = 0 if inlist(h11,0)	
	label variable diarrhea "Diarrhea dichotomous"
	label define diarrhea_label 0 "No diarrhea last 2 weeks" 1 "Diarrhea last 2 weeks"
	label values diarrhea diarrhea_label
	
//	Generate language indicator, defined  if Mayan language spoken at home
	
	gen language = .
	replace language = 1 if inrange(s115,2,23)
	replace language = 0 if inlist(s115,1,96)
	label variable language "Language in home"
	label define language_label 0 "Spanish" 1 "Indigenous language"
	label values language language_label

//	Generate improved sanitation

	gen sanitation = .
	replace sanitation = 1 if inlist(v116,10,11,12,13,14,20,21,22) & inlist(v160,0)
	replace sanitation = 2 if inlist(v116,10,11,12,13,14,20,21,22) & inlist(v160,1)
	replace sanitation = 3 if inlist(v116,23,30,31)
	label variable sanitation "Sanitation"
	label define sanitation_label 1 "Improved sanitation, unshared" 2 "Improved sanitation, shared" 3 "Unimproved sanitation"
	label values sanitation sanitation_label

//	Generate improved drinking water source

	gen water = .
	replace water = 1 if inlist(v113,10,11,12,13,14,51,43)
	replace water = 0 if inlist(v113,20,30,31,32,40,41,42,43,44,61,71,96)
	label variable water "Improved water"
	label define water_label 0 "Not improved water or other" 1 "Improved water"
	label values water water_label
		
//	Birth number

	rename bord birth_order_number

//	Children under 5 in household

	rename v137 children_household
	
	
* _______________ CONFOUNDING VARIABLES ________________________________________

	
// 	Birth spacing variables

	//	Generate indicator variables for antecedent birth spacing

		rename b11 ant_birth_interval_cont
		
	//	Generate indicator variables for subsequent birth spacing

		rename b12 sub_birth_interval_cont
	
//	Confounding variables used in nutritional indicators sensitivity

	//	Dietary diversity 

		//	"alimentos hechos a base de granos, raíces y tubérculos, incluye papilla y alimentos infantiles forti cados de granos"

			gen grains = .
			replace grains = 1 if v412a == 1 | v412b == 1 | v414a == 1 | v414e == 1 | v414f == 1
			replace grains = 0 if v412a == 0 & v412b == 0 & v414a == 0 & v414e == 0 & v414f == 0
			label variable grains "Grains, roots, and tuber"

		//	"legumbres y nueces"

			gen legumes = .
			replace legumes = 1 if v414o == 1 
			replace legumes = 0 if v414o == 0
			label variable legumes "Legumes and nuts"

		//	"fórmula infantil, otras leches diferentes a la materna, queso, yogur u otros productos lácteos"

			gen dairy = .
			replace dairy = 1 if v411 == 1 | v411a == 1 | v414p == 1 | v414v == 1
			replace dairy = 0 if v411 == 0 & v411a == 0 & v414p == 0 & v414v == 0
			label variable dairy "Dairy products"

		//	"carne, pollo, pescado y mariscos (y carnes de vísceras)"
		
			gen flesh = .
			replace flesh = 1 if v414h == 1 | v414m == 1 | v414n == 1
			replace flesh = 0 if v414h == 0 & v414m == 0 & v414n == 0
			label variable flesh "Flesh foods"

		//	"huevos"

			gen eggs = .
			replace eggs = 1 if v414g == 1
			replace eggs = 0 if v414g == 0
			label variable eggs "Eggs"

		//	"frutas y vegetales ricos en vitamina A"

			gen vitamina = .
			replace vitamina = 1 if v414i == 1 | v414j == 1 | v414k == 1
			replace vitamina = 0 if v414i == 0 & v414j == 0 & v414k == 0
			label variable vitamina "Vitamin-A rich foods"

		//	"Otras frutas y vegetales"
			
			gen ofruits = .
			replace ofruits = 1 if v414l == 1
			replace ofruits = 0 if v414l == 0
			label variable ofruits "Other fruits and vegetables"

		//	Calculations

			egen food_groups = rsum(grains legumes dairy flesh eggs vitamina ofruits)
			label variable food_groups "Food groups consumed"
			svy: tab food_groups if inrange(age_child_months,6,23), percent ci
			
			gen diet_diversity = .
			replace diet_diversity = 0 if inrange(food_groups,0,3) & inrange(age_child_months,6,23)
			replace diet_diversity = 1 if inrange(food_groups,4,7) & inrange(age_child_months,6,23)
			svy: tab diet_diversity if inrange(age_child_months,6,23), percent ci
	
	//	Minimum meal frequency
		
		gen meals = .
		replace meals = m39
		replace meals = . if m39 == 8 | inrange(age_child_months,0,5) | inrange(age_child_months,24,59)
		
		gen breastfeeding = .
		replace breastfeeding = v404
		
		gen diet_frequency = .
		replace diet_frequency = 0 if !missing(meals)
		replace diet_frequency = 1 if (breastfeeding == 0 & meals >=4)
		replace diet_frequency = 1 if (breastfeeding == 1 & inrange(age_child_months,6,8) & meals >=2)
		replace diet_frequency = 1 if (breastfeeding == 1 & inrange(age_child_months,9,23) & meals >=3)
		replace diet_frequency = . if inrange(age_child_months,0,5) | inrange(age_child_months,24,59)
		svy: tab diet_frequency if inrange(age_child_months,6,23), percent ci
	
//	Confounding variables used in maternal height imputation-based sensitivity analysis

	//	Maternal absolute height in cm
			
		gen maternal_height = .
		replace maternal_height = v438/10
		replace maternal_height = . if v438 > 9000
		label variable maternal_height "Maternal height in cm"	
		
	//	If missing maternal height variable
		
		gen missing_maternal_height = 0
		replace missing_maternal_height = 1 if v438 > 9000
		
	//	Child HAZ and stunting differences across missing_maternal_height

		ttest haz, by(missing_maternal_height)
		tabulate haz_stunted missing_maternal_height, chi2 column	
			
	
* _______________ ASSESSING VARIABLE LINEARITY AND GENERATING SPLINES __________
	

//	Assessing linearity of HAZ and continuous variables

	//	Child age
	
		preserve
		set seed 111
		sample 6000, count
		*twoway (scatter haz age_child_months, msize(vtiny) jitter(1.5)) (lfit haz age_child_months), name("haz_by_age_child_months",replace)
		restore
		
		//	splines

			mkspline2 age_child_months_spline = age_child_months, cubic nknots(5) displayknots
			regress haz c.age_child_months_spline*
			*adjustrcspline, noci addplot(scatter haz age_child_months, msize(vtiny) jitter(1.0)) name("haz_by_age_child_months_spline",replace)
	
	//	Mother age
	
		preserve
		set seed 111
		sample 6000, count
		*&twoway (scatter haz mother_age, msize(vtiny) jitter(2.5)) (lfit haz mother_age), name("haz_by_mother_age",replace)
		restore
		
		//	Splines

			mkspline2 mother_age_spline = mother_age, cubic nknots(5) displayknots
			regress haz c.mother_age_spline*
			preserve
			set seed 111
			sample 6000, count
			*adjustrcspline, noci addplot(scatter haz mother_age, msize(vtiny) jitter(2.0)) name("haz_by_mother_age_spline",replace)
			restore
		
	//	Wealth
		
		preserve
		set seed 111
		sample 6000, count
		*twoway (scatter haz wealth_index_cont, msize(vtiny)) (lfit haz wealth_index_cont), name("haz_by_wealth_index_cont",replace)
		restore
		
		//	Spline
		
			mkspline2 wealth_index_cont_spline = wealth_index_cont, cubic nknots(5) displayknots
			regress haz c.wealth_index_cont_spline*
			preserve
			set seed 111
			sample 9000, count
			*adjustrcspline, noci addplot(scatter haz wealth_index_cont, msize(vtiny) jitter(2.0)) name("haz_by_wealth_index_cont",replace)
			restore

	//	Maternal height
		
		preserve
		set seed 111
		sample 6000, count
		*twoway (scatter haz maternal_height, msize(vtiny) jitter(5.0)) (lfit haz maternal_height), name("haz_by_maternal_height",replace)
		restore
	
	//	Birth order
	
		preserve
		set seed 111
		sample 6000, count
		*twoway (scatter haz birth_order_number, msize(vtiny) jitter(7.5)) (lfit haz birth_order_number), name("haz_by_birth_order_number",replace)
		restore
		
	//	Children under 5 in household
	
		preserve
		set seed 111
		sample 2000, count
		*twoway (scatter haz children_household, msize(vtiny) jitter(7.5)) (lfit haz children_household), name("haz_by_children_household",replace)
		restore

		
* _______________ REVIEW MISSING DATA __________________________________________	


//	Dependent variables

	mdesc hw5 haz haz_stunted

//	Independent variables

	mdesc method_current method_current_duration method_prior unmet pregnancy_intention

//	Confounding variables

	mdesc child_sex age_child_months mother_age maternal_literacy region rural_urban ethnicity marital_status diarrhea language sanitation water birth_order_number children_household

//	Covariates used in sensitivity analyses

	mdesc ant_birth_interval_cont
	mdesc maternal_height
	
//	Getting a sense about missing data conditional on HAZ missing

	count if !missing(haz) & missing(diarrhea)
	count if !missing(haz) & missing(water)
	count if !missing(haz) & missing(sanitation)
	count if !missing(haz) & missing(maternal_height)

	
* _______________ CONFOUNDER MACROS ____________________________________________

	
//	Generate macro to permit plug-and-play of covariates to models 

	global covariates "1.child_sex c.age_child_months_spline* i.maternal_education_cat i.paternal_education_cat c.mother_age_spline* c.wealth_index_cont_spline* i.maternal_literacy i.region i.rural_urban 1.ethnicity ib0.marital_status 1.diarrhea 1.language ib3.sanitation 0.water i.birth_order_number i.children_household"
	
//	Macro for birth spacing sensitivity analysis

	global covariates_spacing "c.ant_birth_interval_cont"

//	Macro for feeding indicators sensitivity analysis

	global covariates_diet "i.diet_frequency i.food_groups"

	
* _______________ BASIC STATISTICS (TABLE 1) ___________________________________


//	Dependent variables

	svy: mean haz
	estat sd
	svy: tab haz_stunted, percent ci
	svy: tab haz_sev_stunted, percent ci

//	Independent variables
	
	svy: tab method_current_duration, percent ci
	svy: tab method_prior, percent ci
	svy: tab unmet, percent ci
	svy: tab pregnancy_intention, percent ci
		
//	Covariates

	//	Continuous

		svy: mean age_child_months
		estat sd
		svy: mean mother_age
		estat sd
		svy: mean wealth_index_cont
		_pctile wealth_index_cont [iweight=survey_weight], p(25, 50, 75)
	return list
	
	//	Categorical
	
		svy: tab child_sex, percent ci	
		svy: tab rural_urban, percent ci
		svy: tab maternal_education_cat, percent ci
		svy: tab paternal_education_cat, percent ci	
		svy: tab maternal_literacy, percent ci
		svy: tab marital_status, percent ci
		svy: tab region, percent ci
		svy: tab ethnicity, percent ci
		svy: tab language, percent ci
		svy: tab diarrhea, percent ci
		svy: tab sanitation, percent ci
		svy: tab water, percent ci	
	
	//	Categorical variables presented continuously for readability
		
		svy: mean birth_order_number
		estat sd
		svy: mean children_household
		estat sd


* _______________ BIVARIATE ANALYSIS ___________________________________________
	
	
//	Bivariate analysis

	//	HAZ

		svy linearized : mean haz, over(method_current_duration)
		svy linearized : mean haz, over(method_prior)
		svy linearized : mean haz, over(unmet)
		svy linearized : mean haz, over(pregnancy_intention)
			
	//	Stunted
		
		tabout method_current_duration method_prior unmet pregnancy_intention haz_stunted using table1.tab, replace cells(row ci) f(1 1) clab(Row_% 95%_CI) ptotal(none) svy style(tab) lines(none) percent

	//	Severely stunted
	
		tabout method_current_duration method_prior unmet pregnancy_intention haz_sev_stunted using table2.tab, replace cells(row ci) f(1 1) clab(Row_% 95%_CI) ptotal(none) svy style(tab) lines(none) stats(chi2) percent
	

* _______________ PRIMARY REGRESSIONS __________________________________________
	

//	Model set 1: Linear regression - Contraception choice

	//	Model 1a: Current method
	
		svy linearized : regress haz ib0.method_current $covariates, baselevels
		testparm ib0.method_current
				
	//	Model 1b: Current method + prior method
		
		svy linearized : regress haz ib0.method_current ib0.method_prior $covariates, baselevels
		testparm ib0.method_current
		testparm ib0.method_prior
		
	//	Model 1c: Current method (with duration) + prior method
	
		svy linearized : regress haz ib0.method_current_duration ib0.method_prior $covariates, baselevels
		testparm ib0.method_current_duration
		testparm ib0.method_prior

//	Model set 2: Linear regression - Unmet need

	//	Model 2a: Unmet need
		
		svy linearized : regress haz i.unmet $covariates, baselevels
					
	//	Model 2b: Unmet need + pregnancy intention

		svy linearized : regress haz i.unmet ib1.pregnancy_intention $covariates, baselevels
		testparm ib1.pregnancy_intention
		
//	Model set 3: Poisson regression - Contraception choice
		
		//	Model 3A
		
			svy linearized : poisson haz_stunted ib0.method_current_duration ib0.method_prior $covariates, irr baselevels
			testparm ib0.method_current_duration
			testparm ib0.method_prior
		
		//	Model 3B
		
			svy linearized : poisson haz_sev_stunted ib0.method_current_duration ib0.method_prior $covariates, irr baselevels
			testparm ib0.method_current_duration
			testparm ib0.method_prior
						
		//	Model 3C
			
			svy linearized : poisson haz_stunted i.unmet ib1.pregnancy_intention $covariates, irr baselevels
			testparm ib1.pregnancy_intention
				
		//	Model 3D
		
			svy linearized : poisson haz_sev_stunted i.unmet ib1.pregnancy_intention $covariates, irr baselevels
			testparm ib1.pregnancy_intention
		

* _______________ SENSITIVITY 1: PRECEDING BIRTH INTERVAL ______________________


//	Model set 1: Linear regression - Contraception choice

	//	Model 1a: Current method
	
		svy linearized : regress haz ib0.method_current $covariates $covariates_spacing, baselevels
		testparm ib0.method_current
				
	//	Model 1b: Current method + prior method
		
		svy linearized : regress haz ib0.method_current ib0.method_prior $covariates $covariates_spacing, baselevels
		testparm ib0.method_current
		testparm ib0.method_prior
		
	//	Model 1c: Current method (with duration) + prior method
	
		svy linearized : regress haz ib0.method_current_duration ib0.method_prior $covariates $covariates_spacing, baselevels
		testparm ib0.method_current_duration
		testparm ib0.method_prior

//	Model set 2: Linear regression - Unmet need

	//	Model 2a: Unmet need
		
		svy linearized : regress haz i.unmet $covariates $covariates_spacing, baselevels
					
	//	Model 2b: Unmet need + pregnancy intention

		svy linearized : regress haz i.unmet ib1.pregnancy_intention $covariates $covariates_spacing, baselevels
		testparm ib1.pregnancy_intention
		
//	Model set 3: Poisson regression - Contraception choice
		
		//	Model 3A
		
			svy linearized : poisson haz_stunted ib0.method_current_duration ib0.method_prior $covariates $covariates_spacing, irr baselevels
			testparm ib0.method_current_duration
			testparm ib0.method_prior
		
		//	Model 3B
		
			svy linearized : poisson haz_sev_stunted ib0.method_current_duration ib0.method_prior $covariates $covariates_spacing, irr baselevels
			testparm ib0.method_current_duration
			testparm ib0.method_prior
						
		//	Model 3C
			
			svy linearized : poisson haz_stunted i.unmet ib1.pregnancy_intention $covariates $covariates_spacing, irr baselevels
			testparm ib1.pregnancy_intention
				
		//	Model 3D
		
			svy linearized : poisson haz_sev_stunted i.unmet ib1.pregnancy_intention $covariates $covariates_spacing, irr baselevels
			testparm ib1.pregnancy_intention			
	
	
* _______________ SENSITIVITY 2: MATERNAL HEIGHT WITH LISTWISE DELETION_________
			
                                                                                                                             //	Model set 1: Linear regression - Contraception choice

	//	Model 1a: Current method
	
		svy linearized : regress haz ib0.method_current $covariates c.maternal_height, baselevels
		testparm ib0.method_current
				
	//	Model 1b: Current method + prior method
		
		svy linearized : regress haz ib0.method_current ib0.method_prior $covariates c.maternal_height, baselevels
		testparm ib0.method_current
		testparm ib0.method_prior
		
	//	Model 1c: Current method (with duration) + prior method
	
		svy linearized : regress haz ib0.method_current_duration ib0.method_prior $covariates c.maternal_height, baselevels
		testparm ib0.method_current_duration
		testparm ib0.method_prior

//	Model set 2: Linear regression - Unmet need

	//	Model 2a: Unmet need
		
		svy linearized : regress haz i.unmet $covariates c.maternal_height, baselevels
					
	//	Model 2b: Unmet need + pregnancy intention

		svy linearized : regress haz i.unmet ib1.pregnancy_intention $covariates c.maternal_height, baselevels
		testparm ib1.pregnancy_intention
		
//	Model set 3: Poisson regression - Contraception choice
		
		//	Model 3A
		
			svy linearized : poisson haz_stunted ib0.method_current_duration ib0.method_prior $covariates c.maternal_height, irr baselevels
			testparm ib0.method_current_duration
			testparm ib0.method_prior
		
		//	Model 3B
		
			svy linearized : poisson haz_sev_stunted ib0.method_current_duration ib0.method_prior $covariates c.maternal_height, irr baselevels
			testparm ib0.method_current_duration
			testparm ib0.method_prior
						
		//	Model 3C
			
			svy linearized : poisson haz_stunted i.unmet ib1.pregnancy_intention $covariates c.maternal_height, irr baselevels
			testparm ib1.pregnancy_intention
				
		//	Model 3D
		
			svy linearized : poisson haz_sev_stunted i.unmet ib1.pregnancy_intention $covariates c.maternal_height, irr baselevels
			testparm ib1.pregnancy_intention
			

* _______________ SENSITIVITY 3: FEEDING INDICATORS ____________________________		


//	Model set 1: Linear regression - Contraception choice

	//	Model 1a: Current method
	
		svy linearized : regress haz ib0.method_current $covariates $covariates_diet, baselevels
		testparm ib0.method_current
				
	//	Model 1b: Current method + prior method
		
		svy linearized : regress haz ib0.method_current ib0.method_prior $covariates $covariates_diet, baselevels
		testparm ib0.method_current
		testparm ib0.method_prior
		
	//	Model 1c: Current method (with duration) + prior method
	
		svy linearized : regress haz ib0.method_current_duration ib0.method_prior $covariates $covariates_diet, baselevels
		testparm ib0.method_current_duration
		testparm ib0.method_prior

//	Model set 2: Linear regression - Unmet need

	//	Model 2a: Unmet need
		
		svy linearized : regress haz i.unmet $covariates $covariates_diet, baselevels
					
	//	Model 2b: Unmet need + pregnancy intention

		svy linearized : regress haz i.unmet ib1.pregnancy_intention $covariates $covariates_diet, baselevels
		testparm ib1.pregnancy_intention
		
//	Model set 3: Poisson regression - Contraception choice
		
		//	Model 3A
		
			svy linearized : poisson haz_stunted ib0.method_current_duration ib0.method_prior $covariates $covariates_diet, irr baselevels
			testparm ib0.method_current_duration
			testparm ib0.method_prior
		
		//	Model 3B -> this model doesn't run
		
			*svy linearized : poisson haz_sev_stunted ib0.method_current_duration ib0.method_prior $covariates $covariates_diet, irr baselevels
			*testparm ib0.method_current_duration
			*testparm ib0.method_prior
						
		//	Model 3C
			
			svy linearized : poisson haz_stunted i.unmet ib1.pregnancy_intention $covariates $covariates_diet, irr baselevels
			testparm ib1.pregnancy_intention
				
		//	Model 3D -> this model also doesn't run
		
			*svy linearized : poisson haz_sev_stunted i.unmet ib1.pregnancy_intention $covariates $covariates_diet, irr baselevels
			*testparm ib1.pregnancy_intention

				
* ____________________________ SECTION BREAK ___________________________________		
		
		
//	Close log

	log close
